function[] = Func_IterativeParcellation(lhData, rhData, numIter, confidence_threshold, outputDIR, subject, combineLeftRight, Prior_Variability, Prior_SNR)

%lhData should be fs4 surface residual matrix, 2562 x time

% numIter = 6;

% confidence_threshold =2; %this is the ratio between maximum corr and the
% second largest corr

% outputDIR = '/autofs/cluster/nexus/5/users/tianyi/ReDo/'

%subject='bjecog_s03_task_reg';

% combineLeftRight  =1 or 0  ; for a network, do you want to combine the
% activity of both hemispheres for further parcellation?



seedDatalh=[];
seedDatarh=[];

n = size(lhData,1);
m = size(rhData,1);

%%%%%%%%%% Think about the priors, Could include variability and SNR
if range(Prior_Variability) >0
    Prior_Variability =0.4+  0.6*(Prior_Variability- min(Prior_Variability))/(max(Prior_Variability) - min(Prior_Variability)); % normalize the range to 0.4 ~1. Therefore the inv will be between 1~2.5.
end

var_lh = Prior_Variability(1:n);
var_rh = Prior_Variability(n+1:end);

varInv_lh = 1./var_lh;
varInv_rh = 1./var_rh;

if range(Prior_SNR) >0
    Prior_SNR = 0.4+  0.6*(Prior_SNR- min(Prior_SNR))/(max(Prior_SNR) - min(Prior_SNR)) ; % normalize the range to 0.4 ~1. Therefore the inv will be between 1~2.5.
end

SNR_lh = Prior_SNR(1:n);
SNR_rh = Prior_SNR(n+1:end);

%%%%%%%%%%%%%%%%%%%


eval(['mkdir ' outputDIR ]);

for cnt = 1:numIter
    cnt
    eval(['mkdir ' outputDIR '/' subject]);
    eval(['mkdir ' outputDIR '/' subject '/Iter_' num2str(cnt)]);




    if cnt==1
%         [vol, M, mr_parms, volsz] = load_mgh(['/autofs/cluster/ccnl/Hesheng/nexus5/Hesheng/IterativeParcellation//GroupParcNoAnat/lh_network_',num2str(1),'_asym_fs4.mgh']); %ROI_par_40_' num2str(i2) '_lh.mgh']);

         [vol, M, mr_parms, volsz] = load_mgh(['/home/liang/Projects/ASD_QC/Parcellation_template/lh_network_',num2str(1),'_asym_fs4.mgh']); %ROI_par_40_' num2str(i2) '_lh.mgh']);
        ventLh =find(vol>0);
        GrpNetlh{1}= ventLh;
        [vol, M, mr_parms, volsz] = load_mgh(['/home/liang/Projects/ASD_QC/Parcellation_template/rh_network_',num2str(1),'_asym_fs4.mgh']);
        ventRh =find(vol>0);
        GrpNetrh{1}= ventRh;
        for i2=1:18  % get the seed waveforms based on Thomas' parcellation, and weight it by inv(Variability)
            [vol, M, mr_parms, volsz] = load_mgh(['/home/liang/Projects/ASD_QC/Parcellation_template/lh_network_',num2str(i2+1),'_asym_fs4.mgh']); %ROI_par_40_' num2str(i2) '_lh.mgh']);
            idx =find(vol>0);
            seedDatalh(i2,:)= varInv_lh(idx)*lhData(idx,:); % weight the group map using the inverse of individual difference
            GrpNetlh{i2+1} = idx;

            [vol, M, mr_parms, volsz] = load_mgh(['/home/liang/Projects/ASD_QC/Parcellation_template/rh_network_',num2str(i2+1),'_asym_fs4.mgh']); %ROI_par_40_' num2str(i2) '_lh.mgh']);
            idx =find(vol>0);
            seedDatarh(i2,:)= varInv_rh(idx)*rhData(idx,:); % weight the group map using the inverse of individual difference
            GrpNetrh{i2+1} = idx;
            %             seedDatarh(i2,:)=mean(rhData(find(vol>0),:));

        end
        GrpSeedDatalh =seedDatalh;
        GrpSeedDatarh =seedDatarh;

    else

        for i2=1:18  % get the seed waveforms based on the last parcellation
            [vol, M, mr_parms, volsz] = load_mgh([outputDIR '/' subject '/Iter_' num2str(cnt-1) '/NetworkConfidence_' num2str(i2+1) '_lh.mgh']);
            idx = find(vol>=confidence_threshold);

            if length(idx)==0
                maxx = max(max(max(vol)));
                idx = find(vol==maxx);
                disp(['HAAAAHAAAHAAA--Threshold TOOOOOO LARGE for network', num2str(i2)])

            end

            seedDatalh(i2,:)=SNR_lh(idx)*lhData(idx,:); % weight the individual signal based on SNR
            [vol, M, mr_parms, volsz] = load_mgh([outputDIR '/' subject '/Iter_' num2str(cnt-1) '/NetworkConfidence_' num2str(i2+1) '_rh.mgh']);
            idx = find(vol>=confidence_threshold);
            if length(idx)==0
                maxx = max(max(max(vol)));
                idx = find(vol==maxx);
                disp(['HAAAAHAAAHAAA--Threshold TOOOOOO LARGE for network', num2str(i2)])
            end
            seedDatarh(i2,:)=SNR_rh(idx)*rhData(idx,:);

            %             seedDatarh(i2,:)=mean(rhData(find(vol>=confidence_threshold),:));
        end
    end

    %%%%%%%%%%%%%%%% Weight in the group seed in each
    %%%%%%%%%%%%%%%% iteration, should throw in individual variability map as weight in the future%%%%%%%%%%%%%%%

    if cnt>1
        seedDatalh = seedDatalh + GrpSeedDatalh/(cnt-1);
        seedDatarh = seedDatarh + GrpSeedDatarh/(cnt-1);
    end

    %%%%%%%%%%%%%%%%%


    % Do you want to combine the same network of left hemi and right hemi?,
    % If not, uncomment the following 3 lines

    % if combine the two hemisphere, the impact of other hemisphere is decreasing during iteration
    if (combineLeftRight)
        tmp = seedDatalh;
        seedDatalh = seedDatalh+seedDatarh/(cnt+2);
        seedDatarh = seedDatarh+tmp/(cnt+2);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end



    % compute vertex to seed correlation for all vertices

    cValuelh=zeros(size(lhData,1),size(seedDatalh,1));
    cValuerh=zeros(size(rhData,1),size(seedDatarh,1));
%     for i2=1:size(lhData,1)
%
%         for i3=1:size(seedDatalh,1)
%             [R P]=corrcoef(lhData(i2,:),seedDatalh(i3,:));
%             cValuelh(i2,i3)=R(1,2);
%
%
%             [R P]=corrcoef(rhData(i2,:),seedDatarh(i3,:));
%             cValuerh(i2,i3)=R(1,2);
%         end
%     end

    data = [seedDatalh;lhData];
    tmp = corrcoef(data');
    cValuelh = tmp(1:size(seedDatalh,1),end-2561:end)'; %%%2562*seeds

    data = [seedDatarh;rhData];
    tmp = corrcoef(data');
    cValuerh = tmp(1:size(seedDatarh,1),end-2561:end)'; %%%2562*seeds

    cValuelh(isnan(cValuelh(:)))=0;
    cValuerh(isnan(cValuerh(:)))=0;

    %save([outputDIR '/' subject, '/numIter_', num2str(cnt), '/cValue.mat'],'cValuelh','cValuerh');

    %%%%%%%%%%%Further weight in the group map * inv(Variability) by adding
    %%%%%%%%%%%correlation coefficient of 0~ 0.5 according to inv(Variability).

    for i =1:18
        idx = GrpNetlh{i+1};
        cValuelh(idx, i) = cValuelh(idx, i) + (((varInv_lh(idx)-1)/3)/cnt)';

        idx = GrpNetrh{i+1};
        cValuerh(idx, i) = cValuerh(idx, i) + (((varInv_rh(idx)-1)/3)/cnt)';
    end


    %%%%%%%%%%%%%%%%%%%%%%Determine the network membership of each vertex
    %%%%%%%%%%%%%%%%%%%%%%-- Left hemisphere
    %%%%%%%%%%%%%%%%%%%%%%
    data=cValuelh(:,1:18);
    for v=1:size(data,1)

        [cor idx] = sort(data(v,:),'descend');
        parc_membership(v) = idx(1);
        parc_confidence(v) = cor(1)/cor(2);

    end


    for n =1:18
        network = 0*parc_membership;
        confid =0*network;
        network(find(parc_membership==n))=1;
        confid(find(parc_membership==n))=parc_confidence(find(parc_membership==n));

        network(ventLh) = 0; % mask out the ventrical and useless areas in the midline

        save_mgh(network,[outputDIR '/' subject '/Iter_' num2str(cnt) '/Network_' num2str(n+1) '_lh.mgh'],eye(4));
        save_mgh(network.*confid,[outputDIR '/' subject '/Iter_' num2str(cnt) '/NetworkConfidence_' num2str(n+1) '_lh.mgh'],eye(4));
    end

    %%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%Determine the network membership of each vertex
    %%%%%%%%%%%%%%%%%%%%%%-- Right hemisphere

    data=cValuerh(:,1:18);
    for v=1:size(data,1)

        [cor idx] = sort(data(v,:),'descend');
        parc_membership(v) = idx(1);
        parc_confidence(v) = cor(1)/cor(2);


    end


    for n =1:18
        network = 0*parc_membership;
        confid =0*network;
        network(find(parc_membership==n))=1;
        confid(find(parc_membership==n))=parc_confidence(find(parc_membership==n));
        network(ventRh) = 0; % mask out the ventrical and useless areas in the midline

        save_mgh(network,[outputDIR '/' subject '/Iter_' num2str(cnt) '/Network_' num2str(n+1) '_rh.mgh'],eye(4));
        save_mgh(network.*confid,[outputDIR '/' subject '/Iter_' num2str(cnt) '/NetworkConfidence_' num2str(n+1) '_rh.mgh'],eye(4));
    end



    eval(['!cp /home/liang/Projects/ASD_QC/Parcellation_template/lh_network_1_asym_fs4.mgh  ' outputDIR '/' subject '/Iter_' num2str(cnt) '/Network_1_lh.mgh']); %ROI_par_40_' num2str(i2) '_lh.mgh']);

    eval(['!cp /home/liang/Projects/ASD_QC/Parcellation_template/rh_network_1_asym_fs4.mgh  ' outputDIR '/' subject '/Iter_' num2str(cnt) '/Network_1_rh.mgh']); %ROI_par_40_' num2str(i2) '_lh.mgh']);

    eval(['!cp /home/liang/Projects/ASD_QC/Parcellation_template/lh_network_1_asym_fs4.mgh  ' outputDIR '/' subject '/Iter_' num2str(cnt) '/NetworkConfidence_1_lh.mgh']); %ROI_par_40_' num2str(i2) '_lh.mgh']);

    eval(['!cp /home/liang/Projects/ASD_QC/Parcellation_template/rh_network_1_asym_fs4.mgh  ' outputDIR '/' subject '/Iter_' num2str(cnt) '/NetworkConfidence_1_rh.mgh']); %ROI_par_40_' num2str(i2) '_lh.mgh']);


end
